Fit models with metabolic index and temperature trends for each species
Max Lindmark
20 March, 2023
Intro
Main analysis: This scripts reads spatially referenced linear trends over time in biomass, metabolic index and temperature (“05_calculate_trends.Rmd”), and fits spatial sdmTMB models to those
# Load libraries, install if needed
library(tidyverse)
#> Warning: package 'tidyr' was built under R version 4.0.5
library(tidylog)
library(RCurl)
library(devtools)
library(sdmTMB)
library(rMR)
library(patchwork)
library(raster)
library(VoCC)
library(RColorBrewer)
#> Warning: package 'RColorBrewer' was built under R version 4.0.5
library(viridis)
library(ggh4x)
library(ggridges)
# Source code for plots
source_url("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/R/functions/map-plot.R")
#> Warning: package 'sf' was built under R version 4.0.5
theme_set(theme_plot())
pal <- brewer.pal(name = "Dark2", n = 8)Read and plot data
# Read from GH?
d <- read_csv("data/clean/trends.csv")
#> Rows: 92322 Columns: 17
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (3): species, life_stage, id
#> dbl (14): X, Y, quarter, temp_slope, biom_slope, phi_slope, mean_log_biomass...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
glimpse(d)
#> Rows: 92,322
#> Columns: 17
#> $ X <dbl> 320.0167, 320.0167, 320.0167, 320.0167, 320.0167, …
#> $ Y <dbl> 6028.587, 6028.587, 6028.587, 6028.587, 6028.587, …
#> $ species <chr> "cod", "cod", "cod", "cod", "dab", "dab", "dab", "…
#> $ life_stage <chr> "adult", "adult", "juvenile", "juvenile", "adult",…
#> $ quarter <dbl> 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1,…
#> $ temp_slope <dbl> 0.03158410, 0.02985489, 0.03158410, 0.02985489, 0.…
#> $ biom_slope <dbl> -0.5817831, -0.7105018, -0.4165413, -0.4676043, 4.…
#> $ phi_slope <dbl> 0.001582802, 0.004173800, 0.004853177, 0.012797675…
#> $ id <chr> "320.016716551297_6028.58668110674_cod_adult_1", "…
#> $ mean_log_biomass <dbl> 4.3916628, 4.5915365, 2.8446881, 2.9603249, 3.6929…
#> $ mean_temp <dbl> 3.369871, 10.260329, 3.369871, 10.260329, 3.369871…
#> $ mean_oxy <dbl> 8.114926, 6.209476, 8.114926, 6.209476, 8.114926, …
#> $ mean_phi <dbl> 1.5023498, 0.8518787, 4.6064940, 2.6120242, 4.4595…
#> $ mean_log_biomass_ct <dbl> 1.544161015, 1.744034754, -0.002813649, 0.11282313…
#> $ mean_temp_ct <dbl> -2.929083, 3.961374, -2.929083, 3.961374, -2.92908…
#> $ mean_oxy_ct <dbl> 2.3680513, 0.4626013, 2.3680513, 0.4626013, 2.3680…
#> $ mean_phi_ct <dbl> -1.6575600, -2.3080311, 1.4465842, -0.5478856, 1.2…
plot_map +
geom_raster(data = d, aes(X*1000, Y*1000, fill = biom_slope)) +
facet_grid(life_stage ~ species) +
#scale_fill_gradient2() +
scale_fill_viridis()
#> Warning: Removed 2488 rows containing missing values (geom_raster).
plot_map +
geom_raster(data = d, aes(X*1000, Y*1000, fill = temp_slope)) +
facet_grid(life_stage ~ species) +
scale_fill_viridis()
#> Warning: Removed 2488 rows containing missing values (geom_raster).
plot_map +
geom_raster(data = d, aes(X*1000, Y*1000, fill = phi_slope)) +
facet_grid(life_stage ~ species) +
scale_fill_viridis()
#> Warning: Removed 2488 rows containing missing values (geom_raster).
d_cod <- filter(d, species == "cod" & life_stage == "adult")
#> filter: removed 73,974 rows (80%), 18,348 rows remaining
p1 <- plot_map +
geom_raster(data = d_cod, aes(X*1000, Y*1000, fill = biom_slope)) +
scale_fill_gradient2(name = "Biotic trend") +
geom_sf(size = 0.3)
p2 <- plot_map +
geom_raster(data = d_cod, aes(X*1000, Y*1000, fill = temp_slope)) +
scale_fill_gradient2(name = "Temperature trend", breaks = c(0.025, 0.075, 0.125)) +
geom_sf(size = 0.3)
p3 <- plot_map +
geom_raster(data = d_cod, aes(X*1000, Y*1000, fill = phi_slope)) +
scale_fill_gradient2(name = "\u03C6 trend") +
geom_sf(size = 0.3)
((p2 + p3) / (p1 + plot_spacer()))
#> Warning: Raster pixels are placed at uneven horizontal intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven vertical intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven horizontal intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven vertical intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven horizontal intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven vertical intervals and will be
#> shifted. Consider using geom_tile() instead.
ggsave("figures/map_trends_cod.pdf", width = 17, height = 17, units = "cm", device = cairo_pdf)
#> Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted. Consider using geom_tile() instead.
#> Raster pixels are placed at uneven vertical intervals and will be shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven horizontal intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven vertical intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven horizontal intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven vertical intervals and will be
#> shifted. Consider using geom_tile() instead.
# Scale biomass trends (temperature, phi)
d <- d %>%
mutate(group = paste(species, life_stage, sep = "_")) %>%
group_by(group) %>%
mutate(temp_slope_sc = (temp_slope - mean(temp_slope)) / sd(temp_slope),
phi_slope_sc = (phi_slope - mean(phi_slope)) / sd(phi_slope)) %>%
ungroup() %>%
mutate(quarter_f = as.factor(quarter))
#> mutate: new variable 'group' (character) with 8 unique values and 0% NA
#> group_by: one grouping variable (group)
#> mutate (grouped): new variable 'temp_slope_sc' (double) with 92,310 unique values and 0% NA
#> new variable 'phi_slope_sc' (double) with 92,315 unique values and 0% NA
#> ungroup: no grouping variables
#> mutate: new variable 'quarter_f' (factor) with 2 unique values and 0% NAFit sdmTMB models of biomass trends as a function of phi trends
# Make a for loop here..
coefs <- list()
pred_trend1 <- list()
pred_trend2 <- list()
for(i in unique(d$group)){
dd <- filter(d, group == i)
mesh <- make_mesh(dd, c("X", "Y"), n_knots = 200)
plot(mesh)
m <- sdmTMB(
biom_slope ~ phi_slope_sc * mean_phi_ct + mean_log_biomass_ct + quarter_f,
mesh = mesh,
data = dd,
family = gaussian(),
spatial = "on",
spatiotemporal = "off",
control = sdmTMBcontrol(newton_loops = 1)
)
sanity(m)
# Save coefficients??
coefs[[i]] <- tidy(m, conf.int = TRUE) %>% mutate(group = i)
# Now predict the biotic trend for high and low temperature with a sd change in the climate trend
low_mean_clim <- as.numeric(quantile(dd$mean_phi_ct, prob = 0.05))
high_mean_clim <- as.numeric(quantile(dd$mean_phi_ct, prob = 0.95))
mean_log_biomass_ct <- mean(dd$mean_log_biomass_ct)
# Predict
nd <- data.frame(phi_slope_sc = c(-1, 0, -1, 0),
mean_phi_ct = rep(c(low_mean_clim, high_mean_clim), each = 2),
mean_log_biomass_ct = rep(mean_log_biomass_ct, 4),
X = mean(dd$X),
Y = mean(dd$Y),
quarter_f = as.factor(4))
nsim = 1000
pred <- predict(m, newdata = nd, se_fit = TRUE, nsim = nsim)
pred <- data.frame(
mean_phi_ct = rep(c(low_mean_clim, high_mean_clim), each = nsim),
est_change = c(pred[1, ], pred[3, ]),
est_0 = c(pred[2, ], pred[4, ]))
pred <- pred %>%
mutate(bio_trend_delta = est_change - est_0, # if bio trend is positive, bio increases with declines in phi
mean_clim = ifelse(mean_phi_ct == low_mean_clim,
"low_phi", "high_phi"),
group = i)
pred_trend1[[i]] <- pred
# Now make a prediction across a range of trends
# Again I'm predicting for the range of the specific model
nd3 <- data.frame(expand.grid(phi_slope_sc = seq(quantile(dd$phi_slope_sc, probs = 0.05),
quantile(dd$phi_slope_sc, probs = 0.95),
length.out = 50),
mean_phi_ct = c(low_mean_clim, high_mean_clim))) %>%
mutate(mean_log_biomass_ct = mean_log_biomass_ct,
X = mean(dd$X),
Y = mean(dd$Y),
quarter_f = as.factor(4))
pred3 <- predict(m, newdata = nd3, se_fit = TRUE)
pred3 <- pred3 %>%
mutate(lwr = est - 1.96*est_se,
upr = est + 1.96*est_se,
group = i)
pred_trend2[[i]] <- pred3
}
#> filter: removed 73,974 rows (80%), 18,348 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 76,018 rows (82%), 16,304 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 90,406 rows (98%), 1,916 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 90,734 rows (98%), 1,588 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 67,784 rows (73%), 24,538 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 83,990 rows (91%), 8,332 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 83,686 rows (91%), 8,636 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 79,662 rows (86%), 12,660 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
phi_coefs <- dplyr::bind_rows(coefs) %>%
mutate(sig = "N",
sig = ifelse(estimate < 0 & conf.high < 0, "Y", sig),
sig = ifelse(estimate > 0 & conf.low > 0, "Y", sig))
#> mutate: new variable 'sig' (character) with 2 unique values and 0% NA
phi_pred_delta <- dplyr::bind_rows(pred_trend1) %>%
separate(group, c("species", "life_stage"),
sep = "_", remove = FALSE) %>%
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage))
#> mutate: changed 16,000 values (100%) of 'species' (0 new NA)
#> changed 16,000 values (100%) of 'life_stage' (0 new NA)
phi_pred <- dplyr::bind_rows(pred_trend2) %>%
separate(group, c("species", "life_stage"),
sep = "_", remove = FALSE) %>%
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage))
#> mutate: changed 800 values (100%) of 'species' (0 new NA)
#> changed 800 values (100%) of 'life_stage' (0 new NA)
write_csv(phi_coefs, "output/phi_trend_coefs.csv")
write_csv(phi_pred_delta, "output/pred_phi_trend_delta.csv")
write_csv(phi_pred, "output/pred_phi_trend.csv")Plot phi trend models
# Plot coefficients
phi_coefs %>%
filter(term %in% c("phi_slope_sc", "phi_slope_sc:mean_phi_ct")) %>%
ggplot(aes(estimate, group)) +
geom_point() +
geom_vline(xintercept = 0, linetype = 2, color = "grey30") +
facet_wrap(~term, ncol = 1)
#> filter: removed 32 rows (67%), 16 rows remaining
phi_pred_delta2 <- phi_pred_delta %>%
# pivot_wider(names_from = mean_clim, values_from = bio_trend_delta) %>%
# mutate(diff = high_phi - low_phi) %>%
# mutate(Expected = ifelse(diff < 0, "No", "Yes")) %>%
# We expect in general the biotic trends to be negative
# We also expect biotic trends to be MORE negative if it's already in a low phi area
#pivot_longer(c(low_phi, high_phi), names_to = "mean_clim", values_to = "bio_trend_delta") %>%
mutate(mean_clim = ifelse(mean_clim == "high_phi", "High \u03C6", "Low \u03C6")) %>%
mutate(id2 = paste(species, life_stage))
#> mutate: changed 16,000 values (100%) of 'mean_clim' (0 new NA)
#> mutate: new variable 'id2' (character) with 8 unique values and 0% NA
order <- phi_pred_delta2 %>%
group_by(id2, mean_clim) %>%
summarise(mean_delta = mean(bio_trend_delta)) %>%
pivot_wider(names_from = mean_clim, values_from = mean_delta) %>%
mutate(diff = `Low φ` - `High φ`) %>%
arrange(desc(diff)) %>%
mutate(Expected = ifelse(`Low φ` < `High φ`, "Yes", "No"))
#> group_by: 2 grouping variables (id2, mean_clim)
#> summarise: now 16 rows and 3 columns, one group variable remaining (id2)
#> pivot_wider: reorganized (mean_clim, mean_delta) into (High φ, Low φ) [was 16x3, now 8x3]
#> mutate (grouped): new variable 'diff' (double) with 8 unique values and 0% NA
#> mutate (grouped): new variable 'Expected' (character) with 2 unique values and 0% NA
order
#> # A tibble: 8 × 5
#> # Groups: id2 [8]
#> id2 `High φ` `Low φ` diff Expected
#> <chr> <dbl> <dbl> <dbl> <chr>
#> 1 Cod Adult -1.08 -0.608 0.471 No
#> 2 Cod Juvenile -0.248 0.0758 0.323 No
#> 3 Plaice Juvenile -0.0106 0.00798 0.0185 No
#> 4 Flounder Juvenile 0.0106 0.00898 -0.00158 Yes
#> 5 Dab Juvenile 0.0366 -0.0206 -0.0572 Yes
#> 6 Flounder Adult 0.233 0.0883 -0.144 Yes
#> 7 Plaice Adult 0.0617 -0.204 -0.266 Yes
#> 8 Dab Adult 0.404 -0.0714 -0.475 Yes
phi_pred_delta2 <- phi_pred_delta2 %>%
left_join(order %>% dplyr::select(Expected), by = "id2")
#> Adding missing grouping variables: `id2`
#> left_join: added one column (Expected)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 16,000
#> > ========
#> > rows total 16,000
# phi_pred_delta2 %>%
# ggplot(aes(bio_trend_delta, fct_reorder(id2, diff), color = mean_clim, alpha = Expected)) +
# geom_vline(xintercept = 0, linetype = 2, color = "grey30") +
# geom_point(size = 3, position = position_dodge(width = 0.1)) +
# scale_color_manual(values = c("steelblue2", "tomato3")) +
# scale_alpha_manual(values = c(0.4, 1)) +
# labs(y = "", x = "\u0394 biotic trend with a 1 st.dev\ndecline in \u03C6 trend", color = "Mean climate") +
# theme_plot(base_size = 16) +
# theme(legend.direction = "vertical") +
# NULL
phi_pred_delta2$id2 <- factor(phi_pred_delta2$id2, levels = c(order$id2))
phi_pred_delta2 %>%
ggplot(aes(bio_trend_delta, id2, fill = mean_clim, alpha = Expected)) +
geom_vline(xintercept = 0, linetype = 2, color = "grey30") +
geom_density_ridges(color = NA) +
scale_fill_manual(values = c("steelblue2", "tomato3")) +
scale_alpha_manual(values = c(0.4, 0.8)) +
labs(y = "", x = "\u0394 biotic trend with a 1 st.dev\ndecline in \u03C6 trend", fill = "Mean climate") +
theme_plot(base_size = 16) +
theme(legend.direction = "vertical") +
NULL
#> Picking joint bandwidth of 0.00458
ggsave("figures/mi_trend_delta.pdf", width = 17, height = 17, units = "cm", device = cairo_pdf)
#> Picking joint bandwidth of 0.00458Fit sdmTMB models of biomass trends as a function of temperature trends, given mean temperature
# Make a for loop here..
pred_trend_temp <- list()
pred_trend_temp_og <- list()
for(i in unique(d$group)){
dd <- filter(d, group == i)
mesh <- make_mesh(dd, c("X", "Y"), n_knots = 200)
plot(mesh)
m <- sdmTMB(
biom_slope ~ temp_slope_sc * mean_temp_ct + mean_log_biomass_ct + quarter_f,
mesh = mesh,
data = dd,
family = gaussian(),
spatial = "on",
spatiotemporal = "off",
control = sdmTMBcontrol(newton_loops = 1)
)
sanity(m)
# Predict
low_mean_clim <- as.numeric(quantile(dd$mean_temp_ct, prob = 0.05))
high_mean_clim <- as.numeric(quantile(dd$mean_temp_ct, prob = 0.95))
mean_log_biomass_ct <- mean(dd$mean_log_biomass_ct)
nd <- data.frame(expand.grid(temp_slope_sc = c(0, 1),
mean_temp_ct = seq(low_mean_clim,
high_mean_clim,
length.out = 2))) %>%
mutate(mean_log_biomass_ct = mean_log_biomass_ct,
X = mean(dd$X),
Y = mean(dd$Y),
quarter_f = as.factor(4))
nsim <- 1000
pred_og <- predict(m, newdata = nd)
pred <- predict(m, newdata = nd, nsim = nsim)
pred_og <- pred_og %>%
dplyr::select(temp_slope_sc, mean_temp_ct, est) %>%
pivot_wider(names_from = temp_slope_sc, values_from = est) %>%
mutate(delta_bio_trend = `1` - `0`) %>%
mutate(group = i)
pred <- data.frame(
mean_temp_ct = rep(c(low_mean_clim, high_mean_clim), each = nsim),
est_0 = c(pred[1, ], pred[3, ]),
est_change = c(pred[2, ], pred[4, ])) # note different order here than phi
pred <- pred %>%
mutate(bio_trend_delta = est_change - est_0, # if bio trend is positive, bio increases with warming
mean_clim = ifelse(mean_temp_ct == low_mean_clim,
"low_temperature", "high_temperature"),
group = i)
pred_trend_temp[[i]] <- pred
pred_trend_temp_og[[i]] <- pred_og
}
#> filter: removed 73,974 rows (80%), 18,348 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> pivot_wider: reorganized (temp_slope_sc, est) into (0, 1) [was 4x3, now 2x3]
#> mutate: new variable 'delta_bio_trend' (double) with 2 unique values and 0% NA
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 76,018 rows (82%), 16,304 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> pivot_wider: reorganized (temp_slope_sc, est) into (0, 1) [was 4x3, now 2x3]
#> mutate: new variable 'delta_bio_trend' (double) with 2 unique values and 0% NA
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 90,406 rows (98%), 1,916 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> pivot_wider: reorganized (temp_slope_sc, est) into (0, 1) [was 4x3, now 2x3]
#> mutate: new variable 'delta_bio_trend' (double) with 2 unique values and 0% NA
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 90,734 rows (98%), 1,588 rows remaining
#> Warning: The model may not have converged. Maximum final gradient:
#> 0.0107874954216619.
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✖ `b_j` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> ✖ `b_j` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> ✖ `b_j` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> pivot_wider: reorganized (temp_slope_sc, est) into (0, 1) [was 4x3, now 2x3]
#> mutate: new variable 'delta_bio_trend' (double) with 2 unique values and 0% NA
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 67,784 rows (73%), 24,538 rows remaining
#> Warning: The model may not have converged. Maximum final gradient:
#> 0.0211433545542356.
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✖ `b_j` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> ✖ `b_j` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> ✖ `b_j` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> pivot_wider: reorganized (temp_slope_sc, est) into (0, 1) [was 4x3, now 2x3]
#> mutate: new variable 'delta_bio_trend' (double) with 2 unique values and 0% NA
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 83,990 rows (91%), 8,332 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> pivot_wider: reorganized (temp_slope_sc, est) into (0, 1) [was 4x3, now 2x3]
#> mutate: new variable 'delta_bio_trend' (double) with 2 unique values and 0% NA
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 83,686 rows (91%), 8,636 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> pivot_wider: reorganized (temp_slope_sc, est) into (0, 1) [was 4x3, now 2x3]
#> mutate: new variable 'delta_bio_trend' (double) with 2 unique values and 0% NA
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 79,662 rows (86%), 12,660 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> pivot_wider: reorganized (temp_slope_sc, est) into (0, 1) [was 4x3, now 2x3]
#> mutate: new variable 'delta_bio_trend' (double) with 2 unique values and 0% NA
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
pred_trend_temp <- dplyr::bind_rows(pred_trend_temp) %>%
separate(group, c("species", "life_stage"),
sep = "_", remove = FALSE) %>%
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage))
#> mutate: changed 16,000 values (100%) of 'species' (0 new NA)
#> changed 16,000 values (100%) of 'life_stage' (0 new NA)
pred_trend_temp_og <- dplyr::bind_rows(pred_trend_temp_og) %>%
separate(group, c("species", "life_stage"),
sep = "_", remove = FALSE) %>%
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage))
#> mutate: changed 16 values (100%) of 'species' (0 new NA)
#> changed 16 values (100%) of 'life_stage' (0 new NA)
write_csv(pred_trend_temp, "output/pred_temp_trend.csv")pred_trend_temp_og %>%
group_by(group) %>%
mutate(mean_clim = ifelse(mean_temp_ct == min(mean_temp_ct), "Cold", "Warm")) %>%
ungroup() %>%
dplyr::select(mean_clim, group, species, life_stage, delta_bio_trend) %>%
pivot_wider(names_from = mean_clim, values_from = delta_bio_trend) %>%
mutate(diff = Cold-Warm) %>%
mutate(Expected = ifelse(diff > 0, "Yes", "No")) %>%
pivot_longer(c(Cold, Warm), names_to = "mean_clim", values_to = "delta_bio_trend") %>%
mutate(id2 = paste(species, life_stage)) %>%
ggplot(aes(delta_bio_trend, fct_reorder(id2, diff), color = mean_clim, alpha = Expected)) +
geom_point(size = 3, position = position_dodge(width = 0.1)) +
scale_color_manual(values = c("steelblue2", "tomato3")) +
scale_alpha_manual(values = c(0.4, 0.8)) +
geom_vline(xintercept = 0, linetype = 2, color = "grey30") +
labs(y = "", x = "\u0394 biotic trend with a 1 st.dev\nincrease in temperature trend", color = "Mean temperature") +
theme_plot(base_size = 16) +
theme(legend.direction = "vertical") +
NULL
#> group_by: one grouping variable (group)
#> mutate (grouped): new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> ungroup: no grouping variables
#> pivot_wider: reorganized (mean_clim, delta_bio_trend) into (Cold, Warm) [was 16x5, now 8x5]
#> mutate: new variable 'diff' (double) with 8 unique values and 0% NA
#> mutate: new variable 'Expected' (character) with 2 unique values and 0% NA
#> pivot_longer: reorganized (Cold, Warm) into (mean_clim, delta_bio_trend) [was 8x7, now 16x7]
#> mutate: new variable 'id2' (character) with 8 unique values and 0% NA
temp_pred_delta2 <- pred_trend_temp %>%
mutate(mean_clim = ifelse(mean_clim == "low_temperature", "Cold", "Warm")) %>%
mutate(id2 = paste(species, life_stage))
#> mutate: changed 16,000 values (100%) of 'mean_clim' (0 new NA)
#> mutate: new variable 'id2' (character) with 8 unique values and 0% NA
order <- temp_pred_delta2 %>%
group_by(id2, mean_clim) %>%
summarise(mean_delta = mean(bio_trend_delta)) %>%
pivot_wider(names_from = mean_clim, values_from = mean_delta) %>%
mutate(diff = `Warm` - `Cold`) %>%
arrange(desc(diff)) %>%
mutate(Expected = ifelse(Cold > Warm, "Yes", "No"))
#> group_by: 2 grouping variables (id2, mean_clim)
#> summarise: now 16 rows and 3 columns, one group variable remaining (id2)
#> pivot_wider: reorganized (mean_clim, mean_delta) into (Cold, Warm) [was 16x3, now 8x3]
#> mutate (grouped): new variable 'diff' (double) with 8 unique values and 0% NA
#> mutate (grouped): new variable 'Expected' (character) with 2 unique values and 0% NA
order
#> # A tibble: 8 × 5
#> # Groups: id2 [8]
#> id2 Cold Warm diff Expected
#> <chr> <dbl> <dbl> <dbl> <chr>
#> 1 Plaice Adult -0.481 0.261 0.741 No
#> 2 Cod Adult 0.208 0.429 0.221 No
#> 3 Plaice Juvenile -0.000860 -0.00599 -0.00513 Yes
#> 4 Flounder Juvenile -0.0280 -0.0512 -0.0232 Yes
#> 5 Dab Juvenile 0.0172 -0.0241 -0.0413 Yes
#> 6 Cod Juvenile 0.0292 -0.0543 -0.0835 Yes
#> 7 Dab Adult 0.118 -0.214 -0.332 Yes
#> 8 Flounder Adult 0.393 -0.124 -0.518 Yes
temp_pred_delta2 <- temp_pred_delta2 %>%
left_join(order %>% dplyr::select(Expected), by = "id2")
#> Adding missing grouping variables: `id2`
#> left_join: added one column (Expected)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 16,000
#> > ========
#> > rows total 16,000
temp_pred_delta2$id2 <- factor(temp_pred_delta2$id2, levels = c(order$id2))
temp_pred_delta2 %>%
ggplot(aes(bio_trend_delta, id2, fill = mean_clim, alpha = Expected)) +
geom_vline(xintercept = 0, linetype = 2, color = "grey30") +
geom_density_ridges(color = NA) +
scale_fill_manual(values = c("steelblue2", "tomato3")) +
scale_alpha_manual(values = c(0.4, 1)) +
labs(y = "", x = "\u0394 biotic trend with a 1 st.dev\nincrease in temperature trend", fill = "Mean temperature") +
theme_plot(base_size = 16) +
theme(legend.direction = "vertical") +
NULL
#> Picking joint bandwidth of 0.00385
ggsave("figures/temp_trend_delta.pdf", width = 15, height = 17, units = "cm", device = cairo_pdf)
#> Picking joint bandwidth of 0.00385